How Merkel cells transduce mechanical stimuli: A biophysical model of Merkel cells

Merkel cells combine with Aβ afferents, producing slowly adapting type 1(SA1) responses to mechanical stimuli. However, how Merkel cells transduce mechanical stimuli into neural signals to Aβ afferents is still unclear. Here we develop a biophysical model of Merkel cells for mechanical transduction by incorporating main ingredients such as Ca2+ and K+ voltage-gated channels, Piezo2 channels, internal Ca2+ stores, neurotransmitters release, and cell deformation. We first validate our model with several experiments. Then we reveal that Ca2+ and K+ channels on the plasma membrane shape the depolarization of membrane potentials, further regulating the Ca2+ transients in the cells. We also show that Ca2+ channels on the plasma membrane mainly inspire the Ca2+ transients, while internal Ca2+ stores mainly maintain the Ca2+ transients. Moreover, we show that though Piezo2 channels are rapidly adapting mechanical-sensitive channels, they are sufficient to inspire sustained Ca2+ transients in Merkel cells, which further induce the release of neurotransmitters for tens of seconds. Thus our work provides a model that captures the membrane potentials and Ca2+ transients features of Merkel cells and partly explains how Merkel cells transduce the mechanical stimuli by Piezo2 channels.

Merkel cell-neurite complexes generate sustained action potentials with irregular intervals when a sustained mechanical stimulus is applied to the skin [3,6,7].Recent studies found that Piezo2 channels, which are expressed in Merkel cells and Aβ afferents and only generate a transient current, are necessary for the sustained firing of Aβ afferents [7][8][9].However, continuous current injections are required for Aβ afferents to generate sustained action potentials [3,10].Recent theoretical research also thought that rapidly adapting(RA) mechanical-sensitive(MS) channels like Piezo2 cannot account for the sustained responses of Merkel cell-neurite complexes [11].Therefore, how Merkel cells transduce the mechanical stimuli with Piezo2 channels has not been well revealed.
Recent studies found that Merkel cells express many molecules that mediate synaptic vesicle release [12][13][14].Some studies further checked that Merkel cells release neurotransmitter analogs to Aβ afferents [15,16].Merkel cells also generate unique Ca 2+ transients under mechanical stimuli [6,17].Ca 2+ is an important signal that regulates the release of neurotransmitters [18][19][20].Besides, different from other neural cells, Merkel cells do not generate the Na + -related action potentials but have the Ca 2+ and K + dominant membrane potential behaviors [17,[21][22][23].Finally, Piezo2 channels transport ions into the cell under mechanical stimulus, with a preferential transport of Ca 2+ .Therefore, Piezo2 channels, membrane potentials, Ca 2+ transients, and vesicles release together to form a line in Merkel cells to translate the mechanical stimulus.However, a detailed biophysical model of Merkel cells is still lacking.It is still unclear how these parameters participate in the mechanical transduction of Merkel cells.
To address these questions, we establish a biophysical electrophysiological model of Merkel cells.First, we validate our model with experiments of [21][22][23][24].The membrane potentials of Merkel cells under current stimuli vary with the conductances of plasma membrane Ca 2+ channels.The Ca 2+ transients in Merkel cells under high K + solutions and hypotonic shock last for tens of seconds.Then, we explore the roles of ion channels on the plasma membrane, endoplasmic reticulum(ER), and mitochondria(MT) in membrane potentials and Ca 2+ transients regulations of Merkel cells.The results show that Ca v 1.2 channels mainly contribute to the form of the peak of membrane potentials, while K v 1.4 channels inhibit this peak.BKCa, KDR, and Ca v 2.1 channels mainly regulate the steady membrane potentials.Ca v 1.2 channels cause the increase of Ca 2+ concentration at the initial time, while Ca v 2.1 channels, Ryanodine and IP 3 receptors on the ER increase the duration of Ca 2+ transients.Interestingly, the Ca 2+ channels on MT also lengthen the Ca 2+ transients by satisfying the peak of Ca 2+ transients.
Based on the above results, we further study the mechanical transduction of Merkel cells with Piezo2 channels.The results show that Merkel cells generate a continuous neurotransmitter release under a sustained indentation, in which the duration of neurotransmitter release is positively related to indentation depth.This duration lasts for tens of seconds, corresponding to the firing time of Merkel discs under indentation [6,7,24].These results can partly explain how Piezo2 channels(RA MS channels) inspire Merkel cells to have the sustained neurotransmitter excitation on Aβ afferents for their continuous firing.

Kinetics of ions
Most Merkel cells are globular in shape [25], and some Merkel cells in hair follicles have an irregular shape [24].To simplify, here we assume the Merkel cell is a sphere with a radius of r.Merkel cells contain Na + , K + , Cl − , Ca 2+ , and some micro-molecules A − with negative charges inside the cell.Generally, cells keep in a near electro-neutral condition [26,27].We denote n Na , n K , n Cl , n Ca , and n A as the mole numbers of Na + , K + , Cl − , Ca 2+ , and A − inside the cell, respectively.Then their concentrations are where V is the volume of Merkel cells, V = 4/3πr 3 , V ER and V MT are Endoplasmic reticulum and Mitochondria volume.We also set C Na,out , C K,out , C Cl,out , and C Ca,out as Na + , K + , Cl − , and Ca 2+ concentrations in environmental solutions, respectively.Ions inside the cell are regulated by ion channels embedded in the cell membrane.
K + channels.The molecular profiling of Merkel cells shows that they express three kinds of voltage-gated K + channels: K v 4.2, K v 1.4, and K v 8.1 [12].However, the study of Yamashita shows there are only two kinds of ion currents of K + in Merkel cells, which are similar to currents of K v 4.2 and K v 1.4 [21].Further studies elucidated that K v 8.1 channels do not generate ion currents directly but regulate other K + channels' activities indirectly [28].Therefore, here we only consider K v 4.2 and K v 1.4 channels in Merkel cells.
K v 1.4 channels.K v 1.4 channels are activated rapidly with the increase of cell membrane potential.A part of the currents inactivates rapidly, while the remaining currents have a slower inactivation [21,29].The data for modeling K v 1.4 channels was obtained by fitting the data from [29] where g K v 1:4 is the specific membrane conductance of K v 1.4, V m is the membrane potential, E K is the equilibrium (or Nernst) potential of K + , F is the Faraday constant.For m, h fast , and h slow [29], 2 channels are both activated and inactivated rapidly [21,30].The model data was obtained from [21] In addition to the above two K + channels, there are also Ca 2+ -activated K + channels(BK Ca channels) and Delayed-rectifier K channels(KDR channels) in Merkel cells [12,21,23].These two kinds of channels carry the majority of K + currents in Merkel cells [21,23].
BK Ca channels.BK Ca channels are activated by membrane potential and have no obvious inactivation behavior.As the intracellular Ca 2+ concentration increases, the curve of the channels' open probability versus membrane potential has a right shift [23,31].Here we take a similar channel model from [32], where g BK Ca is the specific membrane conductance of BK Ca , E BK Ca is the Nernst potential, For n, where KDR channels present a slow activation behavior, the parameters are taken from [21] (Fig where g KDR is the specific membrane conductance of KDR, E KDR is the Nernst potential, For n, Ca 2+ channels.Ca 2+ plays an important role in Merkel cells' responses under different stimuli like mechanical stimulation, hypotonic shock, etc [12,22,23].Merkel cells mainly express two kinds of voltage-gated Ca 2+ channels: Ca v 1.2 and Ca v 2.1 [12]. Ca v 1.2 channels.Ca v 1.2 channels are L-type(long-lasting) channels, which exhibit a Ca 2+ -dependent inactivation [33,34].The equations were taken from [32], where g Ca v 1:2 is the specific membrane conductance of Ca v 1.2, E Ca is the Nernst potential, for m, h, and hCa, where K hCa = 1uM.Ca v 2.1 channels.Ca v 2.1 channels are P/Q-type channels [33,35], which have no obvious inactivation behavior.The equations were taken from [32], where g Ca v 2:1 is the specific membrane conductance of Ca v 2.1, E Ca is the Nernst potential, for n, Piezo2 channels.Merkel cells cannot translate mechanical stimuli to neural signals without Piezo2 channels, which carry inward currents under the indentation [6,7].Piezo2 channels are mechanical-sensitive(MS) channels, which can also be opened by the microtubule suction of the membrane [36][37][38].Therefore we assume that Piezo2 channels' open probability is regulated by membrane tension.Given that cells have a cortex which is comprised of the cell membrane and a thin, cross-linked actin network lying beneath the membrane, we assume that Piezo2 channels open probability is regulated by the cortex stress σ [39].Piezo2 channels exhibit a fast activation and a fast inactivation [36,40].Besides, Piezo2 channels have another inactivation with a much longer time scale [38].Here we introduce C, O, In, and h slow to represent the channel closed state, open state, short-time inactivation state, and the state beyond long-time inactivation, respectively.Only channels beyond the long-time inactivation state and open state can transport ions.The data for modeling Piezo2 channels was obtained by fitting the data from [36], where J Piezo2 is specific membrane conductance of Piezo2 channels, E Piezo2 is the Nernst potential.Under loading, Under unloading, Other parameters are the same all the time.
The parameters of Piezo2 channels are listed in Table 1, and the specific descriptions of Piezo2 channels are seen in Parameters estimation in S1 Appendix.The currents through Piezo2 channels in simulation and experiments are shown in Fig 1.
Piezo2 channels are non-selective anion channels [36,40].All Ca 2+ , K + , and Na + could flow across Piezo2 channels.But Ca 2+ has the highest priority [36,40].Therefore, we assumed that the currents across Piezo2 channels are Ca 2+ currents.Besides these channels above, Merkel cells have a passive electrical property.We take Na + , K + , Cl − , and Ca 2+ leak channels into consideration.These leak ions flow can be written as Midpoint cortex stress of C (Pa) 1450 [36] σ f1 Reference cortex stress of C (Pa) 120 [36] σ s2 Midpoint cortex stress of the time constant for C (Pa) 1450 [36] σ f2 Reference cortex stress of the time constant for C (Pa) 40 [36] σ s4 Midpoint cortex stress of the time constant for O (Pa) 1450 [36] σ f4 Reference cortex stress of the time constant for O (Pa) 100 [36] σ s7 Midpoint cortex stress of h slow (Pa) 1450 [41] σ f7 Reference cortex stress of h slow (Pa) 200 [41] σ s8 Midpoint cortex stress of the time constant for h slow (Pa) 1450 [41] σ f8 Reference cortex stress of the time constant for h slow (Pa) 2000 [41] g Piezo2 Specific membrane conductance of Piezo2(mS/cm 2 ) 3 [7] E Piezo2 Nernst potential (mV) 6 [36,40] https://doi.org/10.1371/journal.pcbi.1011720.t001where g i,leak and E i are the specific membrane conductance and the equilibrium (or Nernst) potential of i, E i ¼ À RT z i F ln C i C i;out , R is the gas constant, T is the thermodynamic temperature, z i is the valence of ions, i = Na + , K + , Cl − , Ca 2+ .
Pumps, cotransporters, and exchangers.Similar to most cells, Merkel cells also contain ion pumps, cotransporters, and exchangers to keep the ions' balance inside the cell.Here we take Na + /K + pumps, Na + /K + /Cl − cotransporters, K + /Cl − cotransporters, Na + /Ca 2+ exchangers, and Ca 2+ pumps on the plasma membrane into consideration.
Na + /K + pumps, which transport two K + into the cell and three Na + out of the cell once by consuming energy, keep a high K + concentration and a low Na + concentration in the cell [42,43].The intracellular concentrations of Na + and K + affect the rate of Na + /K + pumps [44,45], where P NaKpump is rate constant.K Na,NaK and K K,NaK are constants, K Na,NaK = 10mM, K K,NaK = 140mM [46].Na + /K + /Cl − cotransporters, which transport one Na + , one K + , and two Cl − in the same direction once, mainly loading Cl − into the cell [47][48][49].Their flow is controlled by three ions concentrations [27], where P NKCC1 is the rate constant.K + /Cl − cotransporters, which transport one K + and one Cl − in the same direction once, mainly extruding Cl − from inside the cell [48].Their flow was regulated by the Cl − and K + concentrations across the membrane [50][51][52], where P KCC2 is the rate constant.Plasma membrane Ca 2+ pumps and Na + /Ca 2+ exchangers both help to eliminate Ca 2+ from inside the cell [53,54].A general model of Ca 2+ pumps was used [55], where P Capump is the rate constant, K Capump is constant, K Capump = 0.3uM.Na + /Ca 2+ exchangers exhibit complex dynamic behavior.They are motivated by membrane potential, Ca 2+ and Na + concentrations.Under steady state, Na + /Ca 2+ exchangers transport Ca 2+ outside the cell.When the cell is stimulated by currents or membrane voltage, the direction of Na + /Ca 2+ exchangers flow will be reversed [56,57].The dynamic equation of Na + /Ca 2+ exchangers can be written as [57], where P Cana is the rate constant, K mNa and K mCa are constants, K mNa = 87.5mM,K mCa = 0.5uM, η = 0.1, ksat = 0.35 [57].

Internal Ca 2+ dynamics
Intracellular Ca 2+ sources also play an important role in Ca 2+ regulation in Merkel cells.The elimination of internal Ca 2+ store greatly reduces the Ca 2+ transients [22,23].Endoplasmic reticulum and Mitochondria are the main internal Ca 2+ stores.There are main three channels on the Endoplasmic reticulum (ER) membrane including Ca 2+ ATPase pump, Inositol 1,4,5-trisphosphate receptor (IP 3 receptor), and Ryanodine receptor.Ca 2+ ATPase pumps consume energy to actively load Ca 2+ into endoplasmic reticulum, and their rate increase with cytoplasmic Ca 2+ concentration [58,59] where P pump,ER is the rate constant, K ERpump is the dissociation constant, and their values are listed in Table 2. Ryanodine receptors and IP 3 receptors involve in Ca 2+ -induced Ca 2+ release in Merkel cells [22,23].The Ca 2+ flux by Ryanodine receptor is defined by the equation, where P RYR is the rate constant, K s,RYR and K f,RYR are constant.
τ m time constant (ms) 10000 [22,23] τ h time constant (ms) 20000 [22,23] S ER Surface of ER (um 2 ) 150 [66] V ER Volume of ER (um 3 ) 100 [66] k preIP 3 rate constant (mol/(cm 3 � ms)) 10 −12 (Tuned) S MT Surface of mitochondria (um 2 ) 150 [66] V MT Volume of mitochondria (um 3 ) 10 [66] https://doi.org/10.1371/journal.pcbi.1011720.t002IP 3 receptors have three sites for the combination of Ca 2+ and IP 3 .One site is for activated Ca 2+ , and one site is for inhibitory Ca 2+ [60].It means that the increase of Ca 2+ concentration from a low level, the combination of Ca 2+ to IP 3 receptors activates the receptor, as the Ca 2+ concentration increase to a very high level, the receptor will be inhibited [61].We introduce m and h to represent the activation and inhibitory role of Ca 2+ on IP 3 receptor, where τ m and τ h are time constants.
The remaining site is for IP 3 , and the increase of IP 3 concentration enhances the currents of IP 3 receptors [61].Thus we assume that IP 3 receptors' open probability equals to where P leak and P IP 3 are rate constants, P leak represents the leak flow of Ca 2+ on ER.K IP 3 is the dissociation constant, C Ca,ER is the concentration of Ca 2+ in ER.IP 3 mobilizes Ca 2+ from intracellular stores through the IP 3 receptors [62], and itself also takes a dynamic regulation in a single cell [63].Generally, the entry of Ca 2+ through Ca 2+ channels on the membrane cause the hydrolyzation of phosphatidylinositol 4,5-bisphosphate (PIP 2 ), which produces IP 3 [64].Then IP 3 will convert to other matters even though the Ca 2+ concentration is still high [63].Therefore, we assume the production rate of IP 3 increases with the Ca 2+ concentration and the conversion rate of IP 3 depends on its concentration, where k IP 3 and k dIP 3 are the rate constants, K IP 3 Ca is the dissociation constant, c preIP 3 is the concentration of IP 3 precursors like IP 2 [64].preIP 3 will be also replenished after consumption, and keep a rather stable state [65].So we assume that preIP 3 has a production rate related to its concentration.Then the change of preIP 3 will be According to the above three Ca 2+ channels, the dynamic equation of Ca 2+ in ER will be There are Mitochondrial uniporter(MCU) and mitochondrial Na + /Ca 2+ exchanger (MNCX) on the mitochondrial membrane.MCU takes up Ca 2+ into Mitochondria while MNCX releases Ca 2+ from mitochondria to the cytoplasm.The Mitochondrial uniporter dynamics is controlled by the cytoplasm Ca 2+ concentration [32], where P MCU is the rate constant, K MCU is the dissociation constant.
The dynamics of mitochondrial Na + /Ca 2+ exchanger are controlled by mitochondrial Ca 2+ concentration [32], where The dynamic equation of n Ca,MT will be where S MT is the mitochondria surface, β MT = 0.3, which represents the buffer role for Ca 2+ of mitochondria membrane [32].According to the above ion channels, the ions change inside the cell can be written as where S ref is the reference surface area of Merkel cell.
The membrane potentials of Merkel cells are regulated by currents across total ion channels.The Na + /K + pumps transport two K + into the cell and three Na + out of the cell, so there is one charge out of the cell once.The cotransporters transport a Na + , a K + , and two Cl − in the same direction once, so the total discharge is zero.The Na + /Ca 2+ exchangers transduce three Na + into the cell and one Ca 2+ out of the cell once, so there is one charge out of the cell once.

Then the dynamic equation of
where C m is the specific membrane capacitance.

Volume change
The cell volume V change is dependent on the flow of water across the cell membrane [39] where S is the cell surface area, S = 4πr 2 , J water is water flux, it is controlled by the hydrostatic pressure difference and the osmotic pressure difference across the membrane [39] where α is a rate constant α = 10 −9 cm � ms −1 � Pa −1 , ΔP is the osmotic pressure difference across the membrane, which can be described as Cells have a cortex which is comprised of the cell membrane and a thin, cross-linked actin network lying beneath the membrane.We take the cortex as an elastic layer, the cortex stress satisfies [39] where K is the cortex elastic modulus, K = 6000Pa, S ref represents the cell surface under no stress, σ a is the active contraction stress produced by the actin network, σ a = −100Pa.For a spherical cell, the relationship between the cortex stress and hydrostatic pressure satisfies [39] DP where h c is the cortex thickness, h c = 0.5μm [39].

Validation of the model
There are three main stimuli to study the properties of Merkel cells, which are current pulses, high K + solutions stimuli, and hypotonic shocks [7,[21][22][23].We consider these three kinds of situations below.Current pulses.The rectangular negative current pulses induce the membrane potential dynamics of Merkel cells(Fig 2A and 2B), and the parameters of ion channels are seen in Table 3.The solid lines are simulation results and the triangle lines are experimental results from [21].The results show that Merkel cells have almost passive responses to negative current pulses(Fig 2B).
However, the membrane properties of Merkel cells under positive current pulses are rather different.As shown in Fig 2C and 2D, when the amplitude of the current pulses is small, the membrane potential dynamics of Merkel cells are still like a passive response.As the currents increase(green and wathet blue line), the membrane potentials rise rapidly to a peak, then decrease and gradually stabilize.Unlike most other sensory cells or afferents, Merkel cells do not generate action potentials [7,21].Some Merkel cells exhibit more complex membrane potential dynamics [24] (Fig 3A ).By changing the specific membrane conductance of Ca 2+ and K + channels and rate constants of some pumps on the membrane, we can get similar results.As shown in Fig 3B, the membrane potentials of Merkel cells rapidly reach a step, then fire.Ca 2+ plays an important role in this membrane potential dynamics.If the Ca 2+ concentration is reduced in the external solutions, the fire of membrane potentials disappears(Fig 3C and 3D).
Merkel cells not only have membrane potential regulation but also show Ca 2+ dynamic behaviors.The common stimuli which elicit Ca 2+ in experiments are high K + solutions and hypotonic shock.
High K + stimulus.High K + solutions are made by reducing Na + and increasing K + concentrations in external solutions [22,23].The experimental results from [23] are shown in Fig 4A .Here we increase C K,out by 130mM, and reduce C Na,out by 130mM to represent the high K + solution.The results show that the high K + solution induces a rapid increase in Ca 2+ concentration.After reaching a peak, Ca 2+ concentration decreases slowly(Fig 4B).This Ca 2+ transient could last for tens of seconds in Merkel cells, which is much longer than the timescale of membrane potential dynamics(Fig 2).This Ca 2+ transient is regulated by both plasma membrane Ca 2+ channels and internal Ca 2+ channels in Merkel cells.The inhibition of ER Ca 2+ stores causes a dramatic drop in fluorescence intensity of Ca 2+ (Fig 4C ).By setting the rate Hypotonic shock.The hypotonic shock also results in Ca 2+ transients in Merkel cells.Merkel cells were cultured in solutions that remove a part of Na + and add the same concentration of mannitols.Then hypotonic shock was induced by removing mannitols in solutions [22,23].The hypotonic shock induces Ca 2+ to enter the cell, and the fluorescence intensity of Ca 2+ increase slowly, after reaching a peak, Ca 2+ concentration goes back to the baseline(Fig 5A) [22].In the simulation, hypotonic shock causes the opening of Piezo2 channels, Ca 2+ concentration increases and reaches the peak, then gradually decreases(Fig 5B).But Ca 2+ transients in the experiments of [22] increase slowly while the Ca 2+ transients in the simulation have an immediate increase.An important reason is that hypotonic shock can inhibit the cytoplasmic substances' mobility [70], which hinders the quick entry of Ca 2+ into the cell and diffusion of Ca 2+ in the cell.However, our model didn't take this inhibitory role of hypotonic shock into account, which is a limitation of our model.

Effects of ion channels properties on membrane potentials and Ca 2+ transients
According to experiments from [21][22][23][24], Merkel cells in different locations of skin, different animals, or different experimental conditions have different expression of ions channels.For example, in experiments of yamashita [21], currents across K v 1.4 and K v 4.2 channels with inactivation property are dominated in total K + currents.While in experiments of piskorowski   3.
https://doi.org/10.1371/journal.pcbi.1011720.g003https://doi.org/10.1371/journal.pcbi.1011720.g004[23], BKCa channels carry 50 * 80% of total K + currents.The currents only have a slight decrease with time.Therefore, it is reasonable that these channels have various conductances in different Merkel cells.The results from different experiments and our simulation results are consistent with the opinion that these ion channels and pumps in Merkel cells are the keys that regulate the dynamics of membrane potentials and Ca 2+ transients.However, how these channels control the membrane potentials and Ca 2+ transients is still unclear.Membrane potentials.First, we study how Ca 2+ channels regulate the membrane potentials of Merkel cells.By changing the specific membrane conductance of Ca v 1.2 channels, we find that, when g Cav1. 2 is zero, the membrane potential rapidly increases initially, then reaches a steady state gradually.There is no sharp peak for membrane potential(Fig 6A yellow line).With the increase of g Cav 1.2, the membrane potential rapidly rises to its peak.After a small drop, the membrane potential gradually stabilizes(Fig 6A purple line).The peak value of membrane potentials increases with g Ca v 1:2 (Fig 6A ).
Compared to Ca v 1.2 channels, Ca v 2.1 channels have little role in the peak form of membrane potentials.When g Ca v 2:1 is small, they almost do not affect the membrane potentials(Fig Next, we study the roles of K + channels on membrane potentials.As g K v 1:4 is zero, the membrane potential has a bigger peak(Fig 7A yellow line).With the increase of g K v 1:4 , the peak value of membrane potentials decreases(Fig 7A).This result indicates that K v 1.4 channels inhibit the peak of membrane potentials caused by Ca 2+ channels.
As g K v 4:2 is zero or small, they have little influence on membrane potentials(Fig 7B yellow line).This also means that K v 1.4 channels mainly regulate the membrane potential peak in normal situations.As g K v 4:2 increases to 20, K v 4.2 channels greatly reduce the resting membrane potential, further suppressing the action potential peak(Fig 7B blue line).
When g BKCa is zero, the membrane potential reaches the normal peak but continues to increase over 100mV(Fig 7C yellow line).As g BKCa increases, the membrane potentials decrease to a steady state after reaching the peak(Fig 7C).The steady value of membrane potentials decreases with the increase of g BKCa (Fig 7C ).
When g KDR is zero, the membrane potentials have a little change(Fig 7D yellow line).This result also indicates that BKCa channels mainly control the steady membrane potentials.As g KDR increases, the steady membrane potentials also reduce(Fig 7D).
The changes of membrane potential with different g K v 4:2 .yellow line: g K v 4:2 ¼ 0, purple line: g K v 4:2 ¼ 0:2, green line: g K v 4:2 ¼ 2, blue line: The changes of membrane potential with different g BKCa .yellow line: g BKCa = 0, purple line: g BKCa = 0.18, green line: g BKCa = 1.8, blue line: g BKCa = 18(mS/cm 2 ).(D) The changes of membrane potential with different g KDR .yellow line: g KDR = 0, purple line: g KDR = 0.01, green line: g KDR = 0.1, blue line: g KDR = 1(mS/cm 2 ).https://doi.org/10.1371/journal.pcbi.1011720.g007 In experiments of [21], the inhibition of K + channels and the enhancement of Ca 2+ channels by the external solution containing Ba 2+ make Merkel cells depolarize at smaller currents injection.This is consistent with our results that with the decrease of conductances of K + channels like K v 1.4, BKCa, and KDR channels, the depolarized membrane potentials are more positive(Fig 7A , 7C and 7D).The membrane potentials also form the peak at smaller currents injection, which is consistent with the increase of g Ca v 1:2 and the decrease of g K v 1:4 cause the peak of membrane potentials in simulation (Figs 6A and 7A).
Finally, we also study the influences of internal Ca 2+ receptors and channels on membrane potentials.However, they almost have no function on membrane potentials(Fig 8).
Together, these results indicate that Ca v 1.2 channels mainly help to form the peak of membrane potentials while K v 1.4 channels directly reduce this peak.The coupling between Ca v 2.1, Ca v 1.2 channels, and Ca 2+ pumps contribute to the oscillation of membrane potentials.BKCa and K KDR channels are mainly to maintain a lower steady membrane potential.
Ca  Together, these results indicate that Ca v 1.2 and Ca v 2.1 channels, Ryanodine and IP 3 receptors could increase the Ca 2+ transients, while BKCa and KDR channels reduce the Ca 2+ transients.Combining the inhibition roles of BKCa and KDR on membrane potentials, it can be speculated that BKCa and KDR channels, by controlling the resting membrane potentials, further regulate the voltage-gated Ca 2+ channels to influence the Ca 2+ transients.K v 1.4 and K v 4.2 channels only influence membrane potentials at the initial time, which have less role in Ca 2+ transients.In turn, the depolarization of membrane potentials facilitates the Ca 2+ channels also regulate membrane potentials.Therefore, the membrane potentials and the Ca 2+ transients in Merkel cells are coupled to each other.

How Merkel cells respond to mechanical stimulus
We have studied how ion channels on Merkel cells shape the cell membrane potentials and Ca 2+ transients.However, Merkel cells are mechanical sensory cells, which could transduce tactile stimuli to SA1 afferents [6,8,24].Piezo2 channels are necessary for this transduction [7,9].The SA1 afferents generate continuous action potentials under a static mechanical displacement.Without Merkel cells, SA1 afferents only generate action potentials at the initial moment of stimulation.The knockdown of Piezo2 channels in the Mekel cells causes similar results.A recent study held the view that Piezo2 channels in Merkel cells are not enough to help SA1 afferents generate continuous action potentials [11].Therefore, we want to know how Piezo2 channels and other ion channels in Merkel cells participate in this process.The model is modified as follows.
The changes of Ca 2+ concentration in the cell with different g BKCa .yellow line: g BKCa = 0, purple line: g BKCa = 0.18, green line: g BKCa = 1.8, blue line: g BKCa = 18(mS/cm 2 ).(D) The changes of Ca 2+ concentration in the cell with different g KDR .yellow line: g KDR = 0, purple line: g KDR = 0.01, green line: g KDR = 0.1, blue line: g KDR = 1 (mS/cm 2 ).https://doi.org/10.1371/journal.pcbi.1011720.g010 Exocytosis and endocytosis.Merkel cells transmit information to downstream afferents by releasing neurotransmitters [15,16].Generally, the synapse release can be divided into three parts.First, cell organelles synthesize neurotransmitters into vesicles.Many vesicles form the vesicles pools in the presynapse [71,72].Second, When Ca 2+ ions enter the cell or cytoplasm Ca 2+ concentration increases, vesicles in the pools fuse to the cell membrane(exocytosis) and release neurotransmitters [15,16].Finally, vesicles fused to the membrane recycle into the cell by endocytosis [73,74].Vesicles in the pools keep a relatively fixed number at rest, and the consumption of vesicles can be rapidly replenished.Then we assume that the vesicle synthesis rate is inversely related to the number of vesicles.The exocytosis rate increases with the cytoplasm Ca 2+ concentration [20], and is positively correlated with the number of vesicles.Thus the dynamic equation of vesicles can be written as where n ve is the number of vesicles, k ve and k exo are the rate constants.n ve,s , n ve,f , c Ca,s , and c Ca,f are constants.Endocytosis is regulated by many complex signals, one of which is membrane tension [75][76][77].Generally, when the membrane tension is great, more energy is required for vesicles to form from the cell membrane.Therefore, the endocytosis rate can take the form of where k endo is rate constant.σ s,ve and σ f,ve are constants.Endocytosis reduces the cell membrane, while exocytosis adds the cell membrane.Thus the reference surface S ref changes as where r ve is the average radius of vesicles.The values of parameters are seen in Table 4.
Cell compression.The Merkel cell is a sphere before indentation, it is compressed to a cylinder of changeable radius r as shown in  Reference cortex stress (Pa) 1000 [78] σ ve,f Reference cortex stress (Pa) 15 [78] https://doi.org/10.1371/journal.pcbi.1011720.t004 The above process happens in 100ms.After that, the concentration of Ca 2+ continue to fall (Fig 13D), but the internal Ca 2+ stores start to work.The Ca 2+ in the ER enters the cell through Ryanodine and IP 3 receptors, keeping a relatively high Ca 2+ concentration, and continuously facilitating exocytosis to release neurotransmitters.Exocyotosis can last for tens of seconds.If we assume that neurotransmitters in every vesicle are close.Then the downstream afferents will receive continuous and stable neurotransmitter stimuli, which could generate continuous action potentials.The time of exocytosis increases with the compression depth, which is consistent with the results that the firing time of SAI Afferents increases with the indentation depth in experiments of [6,7].However, if Piezo2 channels were inhibited, no Ca 2+ currents flow into the cell, and the membrane potential is at rest, Merkel cells will not respond to mechanical stimuli(Fig G in S1 Appendix), which is consistent with the results that the knockout of Piezo2 channel or inhibitory of Piezo2 channels by Cd 2+ result in the loss of sustained action potentials in the SAI afferents [7,24].

Discussion
In this article, we develop a biophysically detailed model of Merkel cells for the reception of Merkel cells.We first validate our model of Merkel cells with previous experiments.We then discuss how these ion channels control the membrane potentials and Ca 2+ transients in Merkel cells.Finally, we study how Merkel cells convert the mechanical stimuli into the release of neurotransmitters with the participation of Piezo2 channels, membrane potentials, and Ca 2+ transients.

The responses of Merkel cells under different stimuli
First, the membrane potentials of Merkel cells are controlled by K + and Ca 2+ channels.As shown in Fig 2, Merkel cells exhibit a nearly passive response under small current stimuli.As the injection currents increase, the membrane potential forms a peak and then has a small decrease to the steady state.How do K + and Ca 2+ contribute to this phenomenon?As the conductances of K V 1.4, BKCa, and KDR decrease, the depolarized membrane potential increases, which is consistent with the results that Merkel cells are easier to depolarize with the inhibition of K + channels [21].Oppositely, the increase of the conductances of Ca 2+ channels increases the peak of membrane potentials, which is also similar to the results that the membrane potentials form a peak with the small currents stimulus when Ca 2+ channels were enhanced [21].Ca 2+ also can contribute to the oscillation of membrane potentials.As shown in Fig 3, The membrane potentials of Merkel cells initially have a rapid depolarization, then gradually enter an oscillation.By reducing the Ca 2+ concentration in the external solution, this oscillation disappears [24].In simulation, by changing the conductance of Ca 2+ pumps, there is also the oscillation of membrane potentials.Though there are different Ca 2+ -relative membrane potentials behaviors in Merkel cells.However, it seems to have no obvious qualitative influence on the mechanical transduction of Merkel cells.The action potentials of downstream SAI afferents with two kinds of membrane potentials are still similar [7,24].There may be other roles for the oscillation of membrane potentials of Merkel cells, which needs further study.On the other hand, the oscillation of membrane potentials of Merkel cells is smaller than the amplitude of action potentials of neural cells.The membrane potentials are still limited on low depolarization state by K + channels.
Secondly, a high K + solution is a general way to stimulate Merkel cells, for the main K + channels expressed in the cell membrane.High K + solutions alter the potentials of K + , which causes the flow of K + into the cell, and the membrane depolarizes.The Ca 2+ channels were opened, and Ca 2+ enter the cell.Ca 2+ transients were inspired.This process both link with membrane Ca 2+ channels and internal Ca 2+ stores.The inhibition of internal Ca 2+ stores greatly reduce the amplitude of Ca 2+ transients(Fig 4), which is consistent with the results in [23].The decrease of the conductances of Ca v 1.2 and Ca v 2.1 channels also reduces the Ca 2+ transients(Fig 9 ), with is consistent with that the inhibition of Ca v 1.2 and Ca v 2.1 reduces the Ca 2+ transients in experiments of [12].
Finally, the hypotonic shock also causes the Ca 2+ transients(Fig 5 ), different from the previous two stimuli, hypotonic shock is more like a mechanical stimulus.Piezo2 channels were opened with the increase of cortex stress for the water absorption of the cell.The membrane potential increases, and Ca 2+ enters the cell.It can be seen that the differences of Ca 2+ transients between experiments [22] and simulation(Fig 5).The Ca 2+ concentration in the simulation almost immediately increases, then gradually rises to the peak.However, in the experiments of [22,23], the Ca 2+ transients at the beginning of stimulation increase slowly, then rise and fall.One of the possible reasons is that hypotonic shock can inhibit the cytoplasmic substances' mobility [70], thus hindering the quick entry of Ca 2+ into the cell and the diffusion of Ca 2+ in the cell.However, this role of hypotonic shock has not been confirmed, which is a limitation of our model.

How Merkel cells with Piezo2 channels help Aβ afferents generate sustained action potentials
According to our model, when Merkel cells are compressed, the opening of Piezo2 channels leads to the depolarization of membrane potentials and the influx of Ca 2+ .The increase of membrane potential further causes the influx of Ca 2+ through Ca v 1.2 and Ca v 2.1.Ca 2+ in the cytoplasm causes the exocytosis of neurotransmitters.Due to the rapid inactivation of Piezo2 channels, the membrane potential goes back to the resting state within hundreds of milliseconds.The Ca 2+ channels are also closed.The concentration of Ca 2+ in the cytoplasm decreases.But the internal Ca 2+ store begins to release Ca 2+ into the cell, which can last for tens of seconds.Exocytosis lasts for the same amount of time until Ca 2+ transients disappear.Therefore, Piezo2 channels excite Merkel cells to produce a sustained neurotransmitter release, which further helps Aβ generate prolonged action potentials.So why does the theoretical study of [11] have such a conclusion that RA MS channels like Piezo2 channels cannot account for the sustained responses of Merkel cell-neurite complexes?Because they treat Merkel cells as ordinary nerve cells.A typical feature of these kinds of cells like afferents is that they generate Na +related action potentials.For typical nerve cells, an action potential travels through axons and arrives at the presynapse, which depolarizes the presynaptic membrane potential.Then voltage-gated Ca 2+ channels on the presynaptic membrane are opened, allowing Ca 2+ to enter the cell and activate neurotransmitters release [83,84].But the Ca 2+ transients in these nerve cells quickly return to the baseline [32,85].The number of vesicles released is positively correlated with the number of arriving action potentials.Therefore, the action potentials of downstream nerve cells are also positively related to the action potentials of presynapse [86][87][88].If the nerve cell is compressed, the opening of Piezo2 channels will only make one or several action potentials [46].Therefore, the misuse of neural properties in Merkel cells leads to the wrong conclusion [11].

Limitations of the model
Merkel cell-neurite complexes generate sustained action potentials with irregular intervals when a sustained mechanical stimulus is applied to the skin [3,6,7,24].This irregular interval of action potentials is also one of the important features of the SA1 response.At least in our results, the neurotransmitter release of one Merkel cell is regular.There are at least two aspects that may lead to this phenomenon.First, not only Ca 2+ , IP 3 , but also many other secondary messengers like ATP and cAMP, play important roles in exocytosis in many neurons [89][90][91][92].These secondary messengers interact with each other, regulating the irregular exocytosis, which further contributes to the irregular intervals.Secondly, an afferent is usually linked with several Merkel cells [8,93,94].It means that even if a single Merkel cell releases vesicles regularly, the combination of several or dozens of Merkel cells with small differences in channel properties may lead to irregularities in the whole action potentials.The research of [93] found that the different locations in the skin, cell size, and so on of Merkel cells may determine this irregular interval.Modeling virtual Merkel cells by different membrane capacitance, and different conductance of ion channels also helps downstream afferents generate irregular action potentials.Therefore, further studies on this irregular action potential are needed.
In this article, we didn't consider the diffusion of Ca 2+ and IP 3 in the cytoplasm due to the rather small size of Merkel cells.However, the diffusion of Ca 2+ may cause the Ca 2+ wave within the cell.This means that the Ca 2+ concentration in the cytoplasm will be nonuniform, which may lead to more complex Ca 2+ behaviors.In previous studies, the Ca 2+ transients in the same cell by two same stimuli can be different [22,23].
Although the Piezo2 channels prefer Ca 2+ , Ca 2+ in external solutions is much smaller than Na + and K + .These two kinds of cation channels also flux through Piezo2 channels.Here we ignore ions other than Ca 2+ in Piezo2 channels.It may have quantitative impacts on our results.
Generally speaking, neurotransmitters in vesicles are adequate.But for Merkel cells, the release of neurotransmitters lasts for tens of seconds, and the supplements of vesicles and neurotransmitters may need to be considered.Upon stimulation, the neurotransmitters released in the synapse and vesicles fused to the membrane will be absorbed back into the cell, and some of the neurotransmitters will also diffuse into adjacent solutions [75,95,96].The time scales of reabsorption of neurotransmitters and vesicles may be out of sync.There is even a "kiss and run" approach to neurotransmitter release, where neurotransmitters are released into the synaptic cleft but the membrane of vesicles goes back to the cell immediately.
Therefore, the recycling of neurotransmitters and vesicles may also influence the sustained response of Merkel cells to stimulation.

Conclusion
Although it has been identified that Merkel cells transduce tactile stimuli to SA1 afferents by Piezo2 channels, it remains unclear how Piezo2 channels with fast-inactivation properties help Merkel cells and SA1 afferents generate a slowly adapting response to mechanical stimuli.Here, we develop a biophysically detailed model for the reception of Merkel cells.
We first validate our model with several experimental results [21][22][23][24].Merkel cells exhibit an almost passive response to negative current and small positive current pulses.As the positive current increases, the membrane potential rises to a peak value and then gradually reaches a steady state.Merkel cells also exhibit Ca 2+ transients lasting tens of seconds under high K + solutions and hypotonic shock.
We then discuss how these ion channels control the membrane potentials and Ca 2+ transients in Merkel cells.Our works show that Ca v 1.2 channels contribute to the formation of the peak of membrane potentials, while K v 1.4 channels reduce this peak.Ca v 2.1, BKCa, and KDR channels mainly maintain the steady membrane potentials, Ca v 2.1 channels increase the steady membrane potentials, and BKCa and KDR channels reduce the steady membrane potentials.Interestingly, the oscillations of membrane potentials require the coupling of Ca 2+ channels and Ca 2+ pumps.Compared to these channels on the membrane, the Ryanodine and IP 3 receptors on ER have little role in membrane potentials.Our works also show that Ca v 1.2 channels increase the peak concentrations of Ca 2+ , while Ca v 2.1 channels mainly increase the duration of high Ca 2+ concentrations.Both BKCa and KDR channels suppress the Ca 2+ transients.It can be speculated that these two channels inhibit the depolarization of membrane potentials, thereby inhibiting the voltage-gated Ca 2+ channels on the membrane.Additionally, both Ryanodine and IP 3 receptors on ER increase the Ca 2+ transients.
Finally, we show that Piezo2 channels and internal Ca 2+ stores are sufficient to activate continuous neurotransmitter release to downstream Aβ afferents.Surprisingly, the membrane potentials seem not necessary for Merkel cells to transduce mechanical stimuli to afferents.The depolarization of membrane potentials caused by Piezo2 channels will rapidly repolarize for the inactivation of Piezo2 channels.Oppositely, Ca 2+ flows into the cell through Piezo2 channels, which activates Ryanodine and IP 3 receptors on ER to keep a continuous high Ca 2+ concentration.This Ca 2+ transients facilitate the release of neurotransmitters, and the duration of exocytosis is corresponding to the time of Aβ afferents firing [6,7,24].
Thus, unlike sensory cells that generate Na + -related action potentials, Merkel cells, through stable membrane potentials and Ca 2+ transients regulation, maintain a relatively stable release of neurotransmitters to mechanical stimuli.
However, knowledge about the neurotransmitter interaction between Merkel cells and Aβ afferents is still lacking [97].In further research, we hope to establish a complete model of the Merkel cell and Aβ afferents, which may provide a better description or understanding of the sensing process of the Merkel cell-neurite complexes.

Fig 2 .
Fig 2. The membrane potential dynamics of one Merkel cell under the rectangular current pulses.(A) Negetive current pulses step in 8.14pA (experimental current step in 10pA).(B) Membrane potentials change under negative current pulses.The solid lines are simulation results, and the triangles are experimental results from [21].(C) Positive current pulses step in 8.14pA, Green line: 9 × 8.14pA, blue line: 14 × 8.14pA.(D) Membrane potentials change under positive current pulses.The solid lines are the simulation results, and the triangles are experimental results from[21].The parameters of ion channels are seen in Table3.https://doi.org/10.1371/journal.pcbi.1011720.g002

PFig 3 .
Fig 3.The membrane potential dynamics of one Merkel cell with different properties under the rectangular current pulses.(A) The changes of membrane potential of experimental results from [24].(B) Membrane potentials change under current pulses in simulation, (C) The changes of membrane potential under the condition of reduced Ca 2+ concentration in external solutions from [24].(D) Membrane potentials change under current pulses in simulation at C Ca,out = 5uM.blue line: 27.14pA, red line: 54.28pA, yellow line: 67.85pA.The parameters of ion channels are seen in Table3.

Fig 4 .
Fig 4. The Ca 2+ transients in Merkel cells under high K + stimulation.(A) The fluorescence intensity change of Ca 2+ under high K + solution [23].(B) The cytoplasm concentration of Ca 2+ change under high K + solution: c K,out increases by 130nM and c Na,out reduces by 130nM in simulation.(C) The fluorescence intensity change of Ca 2+ under high K +solution by adding thapsigargin (TG), which is a specific blocker of Ca 2+ pumps on ER, to deplete intracellular Ca 2+ stores[23].(D) The cytoplasm concentration of Ca 2+ change under high K + solution in simulation by setting P RYR = 0, P leak = 0, P IP 3 ¼ 0, and P pump,ER = 0.The parameters of ion channels are seen in Table A in S1 Appendix.

Fig 5 .Fig 6 . 1 . 2 ,
Fig 5.The Ca 2+ transients in Merkel cells under hypotonic shock.(A) The fluorescence intensity change of Ca 2+ under hypotonic shock by removing 30mM mannitol in external solution [22].(B) The cytoplasm concentration of Ca 2+ changes under hypotonic shock by removing 30mM mannitol in external solution in simulation.The parameters of ion channels are seen in Table A in S1 Appendix.https://doi.org/10.1371/journal.pcbi.1011720.g005

S1
Appendix.Description of cell indentation and additional informations about this work.Fig A. The schematic diagram of Merkel cell under indentation.Fig B. K v 1.4 channel.Fig C. K v 4.2 channel.Fig D. KDR channel.Fig E. Piezo2 channel.Fig F. Merkel cell reaches a balanced state given an initial value at rest.Fig G.The responses of Merkel cells under compression.Fig H.The responses of the Merkel cell under three kinds of stimuli with the same parameters.Fig I.The influences of exocytosis rate on vesicle regulation.

Table 3 . Parameters differences of ion channels in Figs 2 and 3. Parameter Value in Fig 2 Value in Fig 3
2+ transients.First, when g Ca v 1:2 is zero, cytoplasmic Ca 2+ concentration only increases slightly under current stimulation(Fig 9A yellow line).As g Ca v 1:2 increases, cytoplasmic Ca 2+ concentration has a rapid rise and fall (Fig 9A).The peak of c Ca increases with g Ca v 1:2 .When g Ca v 2:1 is zero, they don't affect Ca 2+ transients(Fig 9B).As g Ca v 1:2 increases, the Ca 2+ concentration not only has a bigger peak but also keeps at a high level for a longer time(Fig 9B green line).However, if g Ca v 2:1 increases to 0.002, the Ca 2+ concentration rises over to 50μM

Table A .
Ion channel parameters in high K + solutions and hypotonic shock.Table B. Ion concentrations of external solutions.Table C. Initial values of variables.(PDF)

S1 Data. Code data of the model in this work.
(ZIP)